function [t,R,V,R_maneuver,t_maneuver,delVs,...
          maneuverstart_utc,maneuverfinish_utc,maneuvertype,ManeuverMass,...
          ManeuverDuration]...
          = stationkeepingsimulation_day_interval...
            (sim_day,tol_lat,tol_long,...
             target_long,station_keeping_duration_max_hr,day_interval,mass,isp,T)
% Run_For_HPOP_Use
global const Cnm Snm AuxParam eopdata swdata SOLdata DTCdata APdata PC

SAT_Const
constants
load DE436Coeff.mat
PC = DE436Coeff;

% read Earth gravity field coefficients
Cnm = zeros(181,181);
Snm = zeros(181,181);
fid = fopen('GGM03S.txt','r');
for n=0:180
    for m=0:n
        temp = fscanf(fid,'%d %d %f %f %f %f',[6 1]);
        Cnm(n+1,m+1) = temp(3);
        Snm(n+1,m+1) = temp(4);
    end
end
fclose(fid);

% read Earth orientation parameters
fid = fopen('eop19620101.txt','r');
%  ----------------------------------------------------------------------------------------------------
% |  Date    MJD      x         y       UT1-UTC      LOD       dPsi    dEpsilon     dX        dY    DAT
% |(0h UTC)           "         "          s          s          "        "          "         "     s 
%  ----------------------------------------------------------------------------------------------------
eopdata = fscanf(fid,'%i %d %d %i %f %f %f %f %f %f %f %f %i',[13 inf]);
fclose(fid);

% read space weather data
fid = fopen('sw19571001.txt','r');
%  ---------------------------------------------------------------------------------------------------------------------------------
% |                                                                                             Adj     Adj   Adj   Obs   Obs   Obs 
% | yy mm dd BSRN ND Kp Kp Kp Kp Kp Kp Kp Kp Sum Ap  Ap  Ap  Ap  Ap  Ap  Ap  Ap  Avg Cp C9 ISN F10.7 Q Ctr81 Lst81 F10.7 Ctr81 Lst81
%  ---------------------------------------------------------------------------------------------------------------------------------
swdata = fscanf(fid,'%4i %3d %3d %5i %3i %3i %3i %3i %3i %3i %3i %3i %3i %4i %4i %4i %4i %4i %4i %4i %4i %4i %4i %4f %2i %4i %6f %2i %6f %6f %6f %6f %6f',[33 inf]);
fclose(fid);

% read solar storm indices
fid = fopen('SOLFSMY.txt','r');
%  ------------------------------------------------------------------------
% | YYYY DDD   JulianDay  F10   F81c  S10   S81c  M10   M81c  Y10   Y81c
%  ------------------------------------------------------------------------
SOLdata = fscanf(fid,'%d %d %f %f %f %f %f %f %f %f %f',[11 inf]);
fclose(fid);

% read geomagnetic storm indices
fid = fopen('DTCFILE.txt','r');
%  ------------------------------------------------------------------------
% | DTC YYYY DDD   DTC1 to DTC24
%  ------------------------------------------------------------------------
DTCdata = fscanf(fid,'%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d',[26 inf]);
fclose(fid);

% read geomagnetic storm indices
fid = fopen('SOLRESAP.txt','r');
%  ------------------------------------------------------------------------
% | YYYY DDD  F10 F10B Ap1 to Ap8
%  ------------------------------------------------------------------------
APdata = fscanf(fid,'%d %d %f %f %f %f %f %f %f %f %f',[12 inf]);
fclose(fid);

% model parameters
AuxParam = struct('Mjd_UTC',0,'area_solar',0,'area_drag',0,'mass',0,'Cr',0,...
                  'Cd',0,'n',0,'m',0,'sun',0,'moon',0,'sRad',0,'drag',0,...
                  'planets',0,'SolidEarthTides',0,'OceanTides',0,'Relativity',0);
%% Do not Touch
% epoch state (Envisat)
fid = fopen('InitialCondition_data.txt','r');
tline = fgetl(fid);
year = str2num(tline(1:4));
mon = str2num(tline(6:7));
day = str2num(tline(9:10));
hour = str2num(tline(12:13));
min = str2num(tline(15:16));
sec = str2num(tline(18:19));
UTC = [year,mon,day,hour,min,sec];
Y0 = zeros(1,6);
for ii=1:6
    tline = fgetl(fid);
    Y0(ii) = str2num(tline);
end
tline = fgetl(fid);
AuxParam.area_solar = str2num(tline(49:end));
tline = fgetl(fid);
AuxParam.area_drag = str2num(tline(38:end));
tline = fgetl(fid);
AuxParam.mass = str2num(tline(19:end));
tline = fgetl(fid);
AuxParam.Cr = str2num(tline(5:end));
tline = fgetl(fid);
AuxParam.Cd = str2num(tline(5:end));
fclose(fid);

% epoch
Mjd_UTC = Mjday(year, mon, day, hour, min, sec);

AuxParam.Mjd_UTC = Mjd_UTC;
AuxParam.n       = 70;
AuxParam.m       = 70;
AuxParam.sun     = 1;
AuxParam.moon    = 1;
AuxParam.planets = 0;
AuxParam.sRad    = 0;
AuxParam.drag    = 0;
AuxParam.SolidEarthTides = 0;
AuxParam.OceanTides = 0;
AuxParam.Relativity = 0;
Mjd0   = Mjd_UTC;
Step   = 60;   % [s]
We = (2*pi / 86164.1) * [0;0;1];
%[] rotation speed of Earth
%% Change from here

simulation_day = sim_day;
%[] Total simulation days

target_lat = 0;
%[] Longitude and latitude tolerance.

d_day = day_interval;
%[] How many days of integration per while loop ?

N_Step = d_day*24*60*60/Step;
%[] Change d_day to step/

R_target_ecef = Eci2Ecef(UTC)*Y0(1:3)';
V_target_ecef = Eci2Ecef(UTC)*Y0(4:6)' - cross(We,R_target_ecef);
S_target_ecef = [R_target_ecef;V_target_ecef]/1000;
%[]State of the perfect location in ECEF frame.

MAX_STATIONKEEPING_DURATION_HOUR = station_keeping_duration_max_hr;
%[] Search space of station keeping duration

t_save = [];
R_save = [];
V_save = [];
R_maneuver_save = [];
t_maneuver_save = [];
Number_of_maneuver = 0;
delVs = [];
maneuverstart_utc = [];
maneuverfinish_utc = [];
maneuvertype = {};
ManeuverMass = [];
ManeuverDuration = [];
cexh = isp*9.81;
%[] Create  read/write-able memory

currentTime = 0;
% Current Time 

t0 = 0;
%Initial time from given UTC

while currentTime < simulation_day*24*60*60
    
    % Propagation
    Eph = Ephemeris(Y0, N_Step, Step,t0);
    Eph = Eph';
    t = Eph(1,:);
    R = Eph(2:4,:)/1000;
    V = Eph(5:7,:)/1000;
    n = length(t);
    [latitude,longitude] = ECI_Position2LatLong(t,R,UTC);
    
    % Perform maneuver   
    Mjd_UTC_new = Mjday(year, mon, day, hour, min, sec+t(end));
    UTC_new = [year, mon, day, hour, min, sec+t(end)];
    R_target_eci=transpose(Eci2Ecef(UTC_new))*S_target_ecef(1:3);
    V_target_eci=transpose(Eci2Ecef(UTC_new))*S_target_ecef(4:6) + cross(We,R_target_eci);
    S_target_eci = [R_target_eci;V_target_eci];
    t_crossing = t(end);
    [R_maneuver,V_maneuver,t_maneuver,delV_total] = Run_StationKeeping_Linear([R(:,end);V(:,end)],S_target_eci,S_target_ecef,...
        t_crossing,MAX_STATIONKEEPING_DURATION_HOUR,UTC);
    t_save = [t_save,t(:,2:end),t_maneuver];
    R_save = [R_save,R(:,2:end),R_maneuver];
    V_save = [V_save,V(:,2:end),V_maneuver];
    R_maneuver_save = [R_maneuver_save,R_maneuver];
    t_maneuver_save = [t_maneuver_save,t_maneuver];
    Y0 = [R_save(:,end);V_save(:,end)]*1000;
    t0 = t_save(end);
    currentTime = t0(end);
    Number_of_maneuver = Number_of_maneuver+1;
    delVs = [delVs;delV_total];
    maneuverstart_utc = [maneuverstart_utc;CalendarDate(JulianDate([UTC(1:5),UTC(6)+t_maneuver(1)]))];
    maneuverfinish_utc = [maneuverfinish_utc;CalendarDate(JulianDate([UTC(1:5),UTC(6)+t_maneuver(end)]))];
    string = sprintf('%d day interval',day_interval);
    maneuvertype = [maneuvertype;string];
    ManeuverMass = [ManeuverMass,mass*exp(delV_total*1000/cexh)-mass];
    ManeuverDuration = [ManeuverDuration,mass*cexh/T*(1-exp(-delV_total*1000/cexh))];
end
ManeuverDuration = ManeuverDuration';
ManeuverMass = ManeuverMass';
t = t_save;
R = R_save;
V = V_save;
R_maneuver = R_maneuver_save;
t_maneuver = t_maneuver_save;

end